# Title     : TODO
# Objective : TODO
# Created by: Administrator
# Created on: 2019/7/24

library(ggplot2)
library(ggrepel)
library(dplyr)
library(tidyr)
library(tibble)
library(ropls)
library(pROC)
library(egg)
library(purrr)
library(stringr)
library(optparse)

createWhenNoExist <- function(f){
    ! dir.exists(f) && dir.create(f)
}

option_list <- list(
make_option("--i", default = "AllMet.csv", type = "character", help = "metabolite data file"),
make_option("--g", default = "SampleInfo.csv", type = "character", help = "sample group file"),
make_option("--sc", default = "sample_color.txt", type = "character", help = "sample color file")
)
opt <- parse_args(OptionParser(option_list = option_list))

options(digits = 3)

sampleInfo <- read.csv(opt$g, header = T, stringsAsFactors = F) %>%
    select(c("SampleID", "ClassNote"))

head(sampleInfo)

classNotes <- sampleInfo %>%
.$ClassNote %>%
unique()
cn <- combn(classNotes, 2)
for (i in 1 : ncol(cn)) {
    row <- cn[, i]
    group1Name <- row[1]
    group2Name <- row[2]
    parent <- paste0(group2Name, "_", group1Name)
    createWhenNoExist(parent)
    fileName <- paste0(parent, "/Diff_", group2Name, "_", group1Name, ".csv")
    diffName <- read.csv(fileName, header = T, stringsAsFactors = F, comment.char = "", row.names = 1) %>%
        rownames_to_column("Metabolite") %>%
        filter(IS_Final_Pooled_Sig == 1) %>%
        .$Name
    data <- read.csv(opt$i, header = T) %>%
        select(- c("HMDB", "KEGG", "Class")) %>%
        filter(Name %in% diffName) %>%
        gather("SampleID", "Value", - Name) %>%
        spread(Name, "Value") %>%
        inner_join(sampleInfo, by = c("SampleID")) %>%
        filter(ClassNote %in% c(group1Name, group2Name)) %>%
        column_to_rownames("SampleID") %>%
        mutate(ClassNote = factor(ClassNote, levels = unique(ClassNote))) %>%
        as.data.frame()
    print(head(data))
    lrRs <- glm(formula = ClassNote ~ ., data = data, family = binomial(link = "logit"))
    lrRs
    rocRs <- roc(data$ClassNote ~ lrRs$fitted)
    rocRs$auc

    breaks <- seq(0, 1, by = 0.2)

    auc <- round(rocRs$auc, 3)
    pointData <- coords(rocRs, "best", transpose = FALSE) %>%
        round(3) %>%
        as.data.frame()
    p <- ggroc(rocRs, legacy.axes = T, colour = "blue", size = 0.6) +
        theme_bw(base_size = 8.8, base_family = "Times") +
        theme(axis.text.x = element_text(size = 12, vjust = 0.5),
        axis.text.y = element_text(size = 12), legend.position = 'none',
        axis.title.y = element_text(size = 12), legend.margin = margin(t = 0.3, b = 0.1, unit = 'cm'),
        legend.text = element_text(size = 6), axis.title.x = element_text(size = 12), panel.grid.major = element_line(
        linetype = 3, color = "#DDDDDD"), panel.grid.minor = element_blank(), plot.title = element_text(hjust = 0.5,
        size = 12, face = "bold"), plot.margin = unit(c(0.5, 0.5, 0.5, 0.5), "cm")
        ) +
        ggtitle("") +
        geom_abline(intercept = 0, slope = 1, color = "grey", linetype = 1) +
        annotate("text", x = 0.5 + 0.01 , y = 0.5, label = paste0("AUC:", auc), color = "red", hjust = 0, size = 4,
        family = "Times") +
        geom_point(data = pointData, aes(x = 1 - specificity, y = sensitivity), size = 1) +
        geom_text(data = pointData, aes(x = 1 - specificity, y = sensitivity, label = paste0(threshold,
        "(", specificity, ",", sensitivity, ")")), color = "black", hjust = 0, vjust = 1, size = 4, nudge_x = 0.015,
        nudge_y = - 0.015, family = "Times") +
        scale_x_continuous("1 - Specificity", breaks = breaks) +
        scale_y_continuous("Sensitivity", breaks = breaks)

    pdfFileName <- paste0(parent, "/LR_Model_ROC_", group2Name, "_", group1Name, ".pdf")
    ggsave(limitsize = FALSE,pdfFileName, p, width = 6, height = 6)
}








